{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.2 Programming an eigenvalue solver\n",
    "\n",
    "We solve the generalized eigenvalue problem\n",
    "\n",
    "$$\n",
    "A u = \\lambda M u,\n",
    "$$\n",
    "\n",
    "where $A$ comes from $\\int \\nabla u \\nabla v$, and $M$ from $\\int u v$, on the space $H_0^1$.\n",
    "\n",
    "This tutorial shows how to implement linear algebra algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "import math\n",
    "import scipy.linalg\n",
    "from scipy import random\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We setup a stiffness matrix $A$ and a mass matrix $M$, and declare a preconditioner for $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=4, dirichlet=\".*\")\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx\n",
    "pre = Preconditioner(a, \"multigrid\")\n",
    "\n",
    "m = BilinearForm(fes)\n",
    "m += u*v*dx\n",
    "\n",
    "a.Assemble()\n",
    "m.Assemble()\n",
    "\n",
    "u = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse iteration is\n",
    "\n",
    "$$\n",
    "u_{n+1} = A^{-1} M u_n,\n",
    "$$\n",
    "\n",
    "where the Rayleigh quotient\n",
    "$$\n",
    "\\rho_n = \\frac{\\left \\langle A u_n, u_n\\right \\rangle}{\\left \n",
    "\\langle M u_n, u_n\\right \\rangle}\n",
    "$$\n",
    "\n",
    "converges to the smallest eigenvalue $\\lambda_1$, with rate of convergence $\\lambda_1 / \\lambda_2,$ where $\\lambda_2$ is the next smallest eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preconditioned inverse iteration (PINVIT), see [Knyazef+Neymeyr], replaces $A^{-1}$ by an approximate inverse $C^{-1}$:\n",
    "\n",
    "$$\n",
    "\\rho_n = \\frac{\\left \\langle A u_n, u_n\\right \\rangle}{\\left \\langle M u_n, u_n\\right \\rangle} \\\\\n",
    "w_n = C^{-1} (A u_n - \\rho_n M u_n) \\\\\n",
    "u_{n+1} = u_n + \\alpha w_n\n",
    "$$\n",
    "\n",
    "The optimal step-size $\\alpha$ is found by minimizing the Rayleigh-quotient on a two-dimensional space:\n",
    "\n",
    "$$\n",
    "u_{n+1} = \\operatorname{arg} \\min_{v \\in \\{ u_n, w_n\\}} \\frac{\\left \\langle A v, v\\right \\rangle}{\\left \\langle M v, v\\right \\rangle} \n",
    "$$\n",
    "\n",
    "This minimization problem can be solved by a small eigenvalue problem\n",
    "\n",
    "$$\n",
    "a y = \\lambda m y\n",
    "$$\n",
    "\n",
    "with matrices\n",
    "\n",
    "$$\n",
    "a = \n",
    "\\left( \\begin{array}{cc}\n",
    "\\left \\langle  A u_n, u_n \\right \\rangle &\n",
    "\\left \\langle  A u_n, w_n \\right \\rangle \\\\\n",
    "\\left \\langle  A w_n, u_n \\right \\rangle &\n",
    "\\left \\langle  A w_n, w_n \\right \\rangle \n",
    "\\end{array} \\right), \\quad\n",
    "m = \n",
    "\\left( \\begin{array}{cc}\n",
    "\\left \\langle  M u_n, u_n \\right \\rangle &\n",
    "\\left \\langle  M u_n, w_n \\right \\rangle \\\\\n",
    "\\left \\langle  M w_n, u_n \\right \\rangle &\n",
    "\\left \\langle  M w_n, w_n \\right \\rangle \n",
    "\\end{array} \\right).\n",
    "$$\n",
    "\n",
    "Then, the new iterate is\n",
    "\n",
    "$$\n",
    "u_{n+1} = y_1 u_n + y_2 w_n\n",
    "$$\n",
    "\n",
    "where $y$ is the eigenvector corresponding to the smaller eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation in NGSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create some help vectors. `CreateVector` creates new vectors of the same format as the existing vector, i.e., same dimension, same real/ complex type, same entry-size, and same MPI-parallel distribution if any."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = u.vec.CreateVector()\n",
    "w = u.vec.CreateVector()\n",
    "Mu = u.vec.CreateVector()\n",
    "Au = u.vec.CreateVector()\n",
    "Mw = u.vec.CreateVector()\n",
    "Aw = u.vec.CreateVector()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we pick a random initial vector, which is  zeroed on the Dirichlet boundary.\n",
    "\n",
    "Below, the `FV` method (short for `FlatVector`)  lets us access the abstract vector's  linear memory chunk, which in turn provides a \"numpy\" view of the memory. The projector clears the entries at the Dirichlet boundary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r.FV().NumPy()[:] = random.rand(fes.ndof)\n",
    "u.vec.data = Projector(fes.FreeDofs(), True) * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we run the PINVIT algorithm. Note that the small matrices $a$ and $m$ defined above are called `asmall` and `msmall` below. They are of type `Matrix`, a class provided by  NGSolve for  dense matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(20):\n",
    "    Au.data = a.mat * u.vec\n",
    "    Mu.data = m.mat * u.vec\n",
    "    auu = InnerProduct(Au, u.vec)\n",
    "    muu = InnerProduct(Mu, u.vec)\n",
    "    # Rayleigh quotient\n",
    "    lam = auu/muu\n",
    "    print (lam / (math.pi**2))\n",
    "    # residual\n",
    "    r.data = Au - lam * Mu\n",
    "    w.data = pre.mat * r.data\n",
    "    w.data = 1/Norm(w) * w\n",
    "    Aw.data = a.mat * w\n",
    "    Mw.data = m.mat * w\n",
    "\n",
    "    # setup and solve 2x2 small eigenvalue problem\n",
    "    asmall = Matrix(2,2)\n",
    "    asmall[0,0] = auu\n",
    "    asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)\n",
    "    asmall[1,1] = InnerProduct(Aw, w)\n",
    "    msmall = Matrix(2,2)\n",
    "    msmall[0,0] = muu\n",
    "    msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)\n",
    "    msmall[1,1] = InnerProduct(Mw, w)\n",
    "    # print (\"asmall =\", asmall, \", msmall = \", msmall)\n",
    "    \n",
    "    \n",
    "    eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    # print (eval, evec)\n",
    "    u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w\n",
    "    \n",
    "Draw (u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simultaneous iteration for several eigenvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the steps for extending the above to `num` vectors.\n",
    "\n",
    "Declare a `GridFunction` with multiple components to store several eigenvectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num = 5\n",
    "u = GridFunction(fes, multidim=num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create list of help vectors, and a set of random initial vectors in `u`, with zero boundary conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = u.vec.CreateVector()\n",
    "Av = u.vec.CreateVector()\n",
    "Mv = u.vec.CreateVector()\n",
    "\n",
    "vecs = []\n",
    "for i in range(2*num):\n",
    "    vecs.append (u.vec.CreateVector())\n",
    "\n",
    "for v in u.vecs:\n",
    "    r.FV().NumPy()[:] = random.rand(fes.ndof)\n",
    "    v.data = Projector(fes.FreeDofs(), True) * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute `num` residuals, and solve a small eigenvalue problem on a 2 $\\times$ `num` dimensional space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "asmall = Matrix(2*num, 2*num)\n",
    "msmall = Matrix(2*num, 2*num)\n",
    "lams = num * [1]\n",
    "\n",
    "for i in range(20):\n",
    "    \n",
    "    for j in range(num):\n",
    "        vecs[j].data = u.vecs[j]\n",
    "        r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]\n",
    "        vecs[num+j].data = pre.mat * r\n",
    "\n",
    "    for j in range(2*num):\n",
    "        Av.data = a.mat * vecs[j]\n",
    "        Mv.data = m.mat * vecs[j]\n",
    "        for k in range(2*num):\n",
    "            asmall[j,k] = InnerProduct(Av, vecs[k])\n",
    "            msmall[j,k] = InnerProduct(Mv, vecs[k])\n",
    "\n",
    "    ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    lams[:] = ev[0:num]\n",
    "    print (i, \":\", [lam/math.pi**2 for lam in lams])\n",
    "    \n",
    "    for j in range(num):\n",
    "        u.vecs[j][:] = 0.0\n",
    "        for k in range(2*num):\n",
    "            u.vecs[j].data += float(evec[k,j]) * vecs[k]\n",
    "            \n",
    "Draw (u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *multidim-component* select in the *Visualization* dialog box allows to switch between the components of the multidim-GridFunction."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
